Spatially Distributed Stochastic Systems: 
equation-free and equation-assisted preconditioned computation 
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Spatially distributed problems are often approximately modelled in terms of partial differential 
equations (PDEs) for appropriate coarse-grained quantities (e.g. concentrations). The derivation 
of accurate such PDEs starting from finer scale, atomistic models, and using suitable averaging, is 
often a challenging task; approximate PDEs are typically obtained through mathematical closure 
procedures (e.g. mean- field approximations). In this paper, we show how such approximate macro- 
scopic PDEs can be exploited in constructing preconditioners to accelerate stochastic simulations 
for spatially distributed particle-based process models. We illustrate how such preconditioning can 
improve the convergence of equation-free coarse-grained methods based on coarse timesteppers. Our 
model problem is a stochastic reaction-diffusion model capable of exhibiting Turing instabilities. 

PACS numbers: 05.10.Ln, 87.18. Hf, 82.40.Bj, 87.18.La, 82.20.Wt 



I. INTRODUCTION 

Many mathematical models involving reaction and dif- 
fusion (e.g. catalytic reactions, morphogenesis) are based 
on partial differential equations (PDEs) describing the 
evolution of species concentrations in space and timeii^. 
One advantage of such PDE models is the extensive set 
of existing theoretical and computational tools for their 
analysis and efficient simulation. A disadvantage of such 
continuum-based models in certain chemical and biolog- 
ical contexts is the relatively low number of molecules of 
some of the species involved. This may render mean- field 
type PDE models inaccurate; individual-based stochastic 
models become then more appropriate than continuum 
ones. 

Directly using a stochastic, molecular-based model for 
a spatially distributed pattern-forming system will typi- 
cally be very computationally intensive, ft becomes then 
important to extract useful coarse-grained, macroscopic 
information from the microscopic molecular-based model 
using as few detailed simulations as possible. This is the 
goal of equation-free methods^diSiSiLSiSiiS which were de- 
signed for cases where the exact macroscopic equations 
are unavailable in closed form. 

Deriving accurate macroscopic equations rigorously is 
a challenging task (see discussions in Refsiii*iSii^ for fluid 
dynamics, bacteria or eukaryotic cells, respectively). The 
mathematical assumptions leading to closures may not 
be quantitatively correct over large parameter regimes, 
making the PDE models inaccurate there. Equation- 
free methods circumvent the derivation of accurate PDEs 
by using short bursts of fine-scale simulation to esti- 
mate necessary numerical quantities (residuals, action 
of Jacobians) on demand, rather than through explicit 
closed formulas. Here we will illustrate how even approx- 



imate PDE models can be exploited to accelerate "exact" 
(particle-based) simulations. In the context of equation- 
free methods this can be accomplished naturally by us- 
ing the approximate PDEs to construct preconditioners 
in the iterative numerical linear algebra involved in fixed 
point, stability and bifurcation computations. We call 
this procedure "equation-assisted" computation. 

Our model problem is a stochastic reaction-diffusion 
system capable of exhibiting a Turing instabilityi^. Such 
models can serve as a prototypes of more realistic pat- 
tern formation mechanisms during morphogenesis (see 
e.g. Refsii^*i&). The paper is organized as follows: In Sec- 
tion m we introduce this illustrative stochastic reaction- 
diffusion model: the Schnakenberg^'' system of two chem- 
ical species in one dimension. It can predict pattern 
formation under some conditions and it was also used 
previously in complex models of limb development^^. 
In Section IIIII we start by presenting the mean field 
Schnakenberg PDEs at the limit of large particle num- 
bers, and briefly summarize their bifurcation behavior 
in a parameter regime where they exhibit pattern forma- 
tion. We briefly describe, in Sect ion HV1 the computation 
of bifurcation diagrams using a timestepper based ap- 
proach - both deterministic and "equation-free", based 
on a stochastic simulator implementing a spatially dis- 
cretized version of the Gillespie Stochastic Simulation 
AlgorithmiSi (SSA) . We also discuss basic features of pre- 
conditioning and "equation-assisted" bifurcation compu- 
tations. We then present, in Section our "equation- 
assisted" results and discuss their comparison with the 
"equation-free" case. Here the mean field PDE is used to 
construct a preconditioner, to accelerate the numerical 
linear algebra in our coarse-grained steady state compu- 
tations. We conclude with a brief summary and discus- 
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II. THE STOCHASTIC REACTION-DIFFUSION 
MODEL 

We consider the Schnakenbergl^ system of two chem- 
ical species U and V with the following reaction mecha- 
nism: 



A ^ U ^ 



C 



B 

2U + V 



V 



(1) 

(2) 
(3) 



Here Eq. describes production and degradation of U 
and Eq. ||2Jl describes production of V. Moreover, U is 
also produced in the reaction Eq. We will assume 
that the concentrations of A and B are constants. 

To simulate stochastically Eqs. (QJ - ©, one can use 
the Gillespie SSA, a standard way to model stochasti- 
cally a spatially homogeneous (well mixed) chemical sys- 
tem. The algorithm is based on answering two essential 
questions at each time step: when will the next chem- 
ical reaction occur, and what kind of reaction will it 
be? Gillespiei^ derived a simple way to answer these two 
questions - at each step, the computer performs a reac- 
tion, updates numbers of reactants and products and con- 
tinues with another time step until the algorithm reaches 
a time of interest. 

Next, we introduce diffusion to the system. We assume 
that U diffuses with (macroscopic) diffusion coefficient 
di and V diffuses with (macroscopic) diffusion coefficient 
d2- We consider a spatially one-dimensional domain - 
the interval [0, 1] with suitable boundary conditions as 
specified later. The generalization of Gillespie's ideas to 
spatially nonhomogeneous systems can be found in the 
literature (see e.g. Refsi22*^). Here, we follow the most 
straightforward way, adding diffusion as another set of 
"reactions" to the system. Namely, we divide our domain 
into m boxes (small intervals) of length h = 1/m. We 
denote by Ui and Vi the number of respective molecules 
in the spatial interval [{i — l)/m,i/m] for i = 1, . . . ,to. 
This means that we describe the state of the stochastic 
reaction-diffusion system by two m-dimensional vectors 
U = [C/i, [/2 . . . , [/™] , V = [^1, ^2 . . . , Kn] and we con- 
sider the following reactions at each time step 




A 
B 



2U, + V, 



1 = 1,. 



i = f , . 

1 = 2,. 



,m, 



(4) 



, m-f, (5) 

. ,TO, (6) 



where Eq. Q means that we implement the Schnaken- 
berg reaction mechanism (Eqs. - (O) in every spa- 
tial box. Equations ® - © describe diffusion; the 
transition rates between boxes are denoted by di and 
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FIG. 1: (Color online) Solutions_of Eqs. © - %_w\th pe- 
riodic boundary conditions for di = 5x 10~*, d2 = 0.01, 
0.03 and 0.06 aX t = 100 with an initial condition being the 
perturbed uniform steady state. Spatial patterns develop for 
d2 > 0.0198; the solutions were shifted to have a local U min- 
imum at a; = 0. We plot U (red curve) and V (blue curve) in 
the same picture. 



^2 and they are connected to the macroscopic diffusion 
coefficients through the formulas di = di/h? = dim?, 
and d2 = d^jh^ = d^rr? ■ The Eqs. |@J - together 
with suitable boundary conditions, will be simulated us- 
ing Gillespie SSA and will be our illustrative stochastic 
reaction-diffusion problem in this paper. 



III. DETERMINISTIC ANALYSIS OF THE 
MODEL PROBLEM 

If we have enough molecules in the system, then 
Eqs. Q - © are well approximated by a system of 
two reaction-diffusion PDEs for the species concentra- 
tions U and V; at the mesh points Xi, U{xi) = Ui/u, 
V(x.i) = Vi/uj, where Xi = [i ~ 1/2) ft,. The constant w 
can be interpreted as the number of molecules in the box 
corresponding to a dimensionless concentration of 1. Re- 
action and diffusion rates are scaled as follows: A = luA, 
B = LuB, ki = fci, k2 = k2, ks = k^, k^ = k^jJ^, 
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FIG. 2: (Color online) Comparison of histograms obtained by 
stochastic simulations with the deterministic results given by 
solution of Eqs. Q-® with periodic boundary conditions. 
Solutions were shifted to have the local U minimum at s = 0. 



di — di/K^, d2 = d2/h?^ where A and B are (constant) 
concentrations of reactants A and B, kj, j = 1,2,3,4, 
are macroscopic reaction rate constants and c?i, c?2 are 
macroscopic diffusion coefficients. Next (instead of fur- 
ther nondimensionalization) we simply choose values for 
the six kinetic parameters and we study the Turing insta- 
bilities of the resulting model. In the rest of this paper, 
we put A = 1, B = 1, ki = 1, k2 = 2, ks — 3, k4 — 1. 
Then the scaling of kinetic constants reads as follows 

A^UJ, B^LU, fci = 1, fc2 2, 



1 



di 



fcs ==3, ki^ di = — J, ^2 = i-TT 



(7) 



Passing the number lu of molecules in a box to infinity and 
the box length h to zero, i.e. uj —^ oo and /i — > 0+, one 
can derive the following system of macroscopic partial 
differential equations for concentrations U and V. 



dU — d^U 



dV_ 

dt 



2U + U^V 



-u'v 



(8) 



(9) 



Here, U : [0,1] -> [0,oo), V : [0,1] ^ [0, oo) and suit- 
able boundary conditions (e.g. no-flux, periodic) must 



be introduced. Considering no-ffux or periodic boundary 
conditions, one can easily verify that the homogeneous 
steady state of Eqs. - © is given by Uh{x,t) = 2, 
Vh{x,t) = 3/4. Linearizing Eqs. l|HJl - ®i one sees that 
the homogeneous steady state is stable for c?i = c?2 = 0, 
i.e. when no diffusion is present in the system. In fact, 
the same result holds whenever di = d2'. no spatial pat- 
terning can be expected if the diffusion coefficients of 
both species are the same. However, Turingi^ showed 
that the homogeneous steady state {Uh,Vh) might be- 
come unstable for di 7^ ^2- Indeed, linearizing Eqs. (jHl) 
~ ©, one finds that the steady state {Uh,Vh) is unsta- 
ble and spatial patterns develop if (^2 > 39.6 di. In this 
paper, we fix the diffusion coefficient di = 5 x 10^^. Spa- 
tial patterns may then develop for d2 > 0.0198. We show 
numerically computed solutions of Eqs. ||HJ) - © with 
periodic boundary conditions for different values of the 
diffusion coefficient d2 in Fig.^ The initial condition was 
chosen to be (Uh,Vh) with small additive random noise. 
The graphs of U (red solid curve) and V (blue dashed 
curve) are plotted at dimensionaless time t = 100, and 
can be practically considered as steady states; this has 
been confirmed also through a steady state solver. 

In Fig. 12 we compare representative SSA results with 
the deterministic ones. We divide the domain [0, 1] into 
m = 200 boxes and using lu = 100 (i.e., defining density 
scaling so that dimensionless density 1 corresponds to 
100 molecules of the relevant chemical species in a box) 
and we choose the values of the parameters as in Eq. (0 
together with di = 5 x 10~* and di — 0.06. We will 
quantify the fluctuations of the stochastic simulations at 
stationarity below. 



IV. EQUATION-ASSISTED COMPUTATION: 
THEORETICAL FRAMEWORK 



A. Numerical bifurcation computations 

The computer-assisted study of Turing patterns in a 
deterministic PDE context requires the numerical com- 
putation and parametric continuation of steady states. 
Spatially distributed PDE steady states in a bifurcation 
diagram are, in general, computed by discretizing the 
PDE into a (sufficiently) large set of ODEs of the type 
dx/dt = f(x;p), finding the roots of f(x;p) = and 
continuing them in parameter space. Here x S M.^ is a 
vector containing the system state (the discretized con- 
centrations of U and V), f : M.^ M.^ is the right hand 
side of the discretization of Eqs. © and I© and p e K*^ 
is a M-dimensional parameter vector; here M = 1 since 
we consider the single parameter d2. Pseudo-arclength 
continuation and branch switching are by now standard 
numerical tools that have been incorporated in special 
purpose packages like AUTOS^ or CONTENT^. 

For deterministic problems for which a good dynamic 
simulator is available, the so-called "timestepper-based" 
approach allows the computation of bifurcation diagrams 
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FIG. 3: (Color online) Bifurcation diagram of the determin- 
istic Schnakenberg system of two species {U and V) with re- 
spect to the diffusion coefficient d2. The diffusion coefficient 
for U, di is fixed at 5 x fO~*. Axis notation: a and b are 
the first Fourier coefficients of the solution (for sin(2;) and 
cos(3;) respectively). Different steady state solution branches 
are plotted in different colors. Solid (dashed) lines represent 
stable (unstable) steady states. The bifurcation points along 
the solution curves are marked by circles and denoted with the 
corresponding bifurcation type ("SN" for Saddle-Node bifur- 
cation and "PF" for Pitchfork bifurcation). The straight line 
represents the uniform steady state (as shown in the solution 
profile for point 13). The inset is a blowup of the bifurcation 
diagram for small d2. Representative solution profiles (solid 
(dashed) line for U {V)) on different branches (numbered and 
marked by "x") are also included. 



in the form of a "wrapper" around the dynamic simula- 
tor (see e.g. Reism^^). Given the current state x as 
an initial condition, the timestepper computes the fu- 
ture (after a "reporting time" T) state $t(x; p) = x(r). 
Steady states of the original system are then found as 
fixed points of 4'(x) = where '4'(x) = x — $t(x; p). 

The bifurcation diagram in Fig. |3| has been computed 
in both ways (giving, of course, identical results); a dis- 
cussion of some pertinent details can be found in the 
Appendix. Fixed point algorithms, like the Newton- 
Raphson iteration, constitute the workhorse of these so- 
lutions of large sets of coupled, nonlinear algebraic equa- 



tions; these involve the repeated solutions of large sets of 
linear, coupled algebraic equations. Consider now per- 
forming these repeated linear solves through matrix-free 
iterative linear solvers (such as GMRES2&); in a nonlin- 
ear equation context we will typically use a matrix-free 
Newton-GMRES solver—. For the timestepper-based 
computation (solving the nonlinear system ^(x) = 0) we 
do not need to compute the Jacobian = d^{x)/dx. 
We only need to compute matrix- vector products of this 
Jacobian with given known vectors v, which can be es- 
timated by a finite different approximation "D"^ ■ v « 
[\I'(x -I- ev) — ^'(x)]/e with suitably small e. 

Such matrix-free linear algebra methods constitute an 
important component of equation-free bifurcation calcu- 
lations. In this context macrosocopic, coarse-grained 
equations are not explicitly available; yet we believe they 
exist, and we do have available a fine-scale (in this case, 
stochastic) dynamic simulator. We can then substitute 
the (unavailable) deterministic timestepper with a fine 
scale, stochastic timestepper involving lifting, evolving, 
and restriction steps^iSiSi. 

This provides an estimate of the (unavailable) deter- 
ministic timestepper for the (unavailable) closed macro- 
scopic evolution equations, obtained on demand through 
the stochastic simulator. All computations of the matrix- 
free Newton-GMRES involve calls to such a timestep- 
per (with systematically chosen initial conditions); the 
stochastic simulator can then be used to numerically 
compute coarse-grained bifurcation diagrams such as the 
one in Fig. |21 even in the absence of closed macroscopic 
evolution equations (i.e., equation-free). 

In this paper, we propose an "equation-assisted" ap- 
proach that accelerates equation-free computations by 
linking approximate deterministic models with accurate 
stochastic ones: equation-free bifurcation computations 
(based on the coarse timestepper) are preconditioned us- 
ing the timestepper of an (approximate) deterministic 
model. 



B. Preconditioning and Equation- Assisted 
Computation 

Good discussions of the basic features of Newton- 
GMRES, as well as pseudocode and MATLAB imple- 
mentations can be found in Refs— (2^. Consider solving 
a general set of N nonlinear equations with N unknowns, 
'3i'(x) ~ 0; the linear equations to be solved at each New- 
ton step are of the form 



AAx = b 



for A e R'^'<^, b e M^, x e with 



(10) 



A = X>*U=x., b= *(x,), 

where Xc is the current solution guess at each Newton 
step. For every iterative linear solve, it is important to 
note that, at each iteration in GMRES, only one call to 
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'i'(x) is needed. GMRES does not require the Jacobian 
matrix X>\I'|x=Xc to be computed explicitly. The Jaco- 
bian matrix always occurs in the form of a matrix-vector 
product, which can be approximated by finite differences: 
Avfc w [\I'(xc+/iVfc) — '4'(xc)]//i where h is suitably small. 

When Newton-GMRES is used in steady state com- 
putations using the coarse timestepper (i.e. ^'**(x) = 
X— $|?(x; p) = 0, where $|?(x; p) is the coarse timestep- 
per based on the stochastic simulator), each evaluation of 
*'**(x) involves evolving the coarse timestepper #^*(x; p) 
for time T, which is often computationally intensive; pos- 
sibly several replica simulations need to be performed for 
variance reduction. It thus becomes an important task to 
reduce the total number of function evaluations to con- 
vergence. As discussed in Ref.^^, GMRES requires less 
overall function evaluations when the eigenvalues of the 
matrix (i.e., A in Eq. H10() ') are more clustered. For a 
given linear system in the form of Eq. H1(J|) , the precon- 
ditioning of GMRES involves finding a regular matrix P, 
such that the preconditioned linear system, 

PAAx = Pb (11) 

leads to a more clustered eigenvalue spectrum. Solving 
Eq. H10|l is equivalent to solving Eq. Hll() . It is well known 
that the system in Eq. Hll|l will have better properties 
(from a numerical point of view) than the original sys- 
tem in Eq. H1U|I if P is close to the inverse of A, i.e. if 
||P — A~^|| is small using a suitable matrix norm. Hence, 
preconditioning by an appropriate matrix P can improve 
the efficacy of numerical solvers for Eq. 11U|) : the goal of 
this paper is to show how this preconditioning idea can 
be applied to equation-free stochastic reaction-diffusion 
problems (and spatially distributed evolution problems 
more generally). 

An approximate deterministic evolution equation for 
the stochastic system statistics may be available (Eqs. © 
and lO), based on closure assumptions, which is not ac- 
curate enough to compute with; yet we can take advan- 
tage of such an evolution equation by using it to cre- 
ate a "good" preconditioning matrix P. This precon- 
ditioning is implemented here by multiplying the orig- 
inal output of each (stochastic simulation based) eval- 
uation of *'**(x) with P = (X"4'''''*|x=xJ~\ where 
^dei^^-j ^ X — $^'^*(x; p) is defined using the determin- 
istic mean field PDE timestepper $^'^*(x;p). That is, 
we use the deterministic timestepper of the approximate 
PDEs (here, at the current solution guess Xc) to help ac- 
celerate the equation-free Netwon-GMRES computation; 
a much simpler preconditioning scheme would constantly 
use the inverse of the deterministic timestepper at the de- 
terministic steady state. At each Newton step the linear 
equation set to be actually solved by GMRES after pre- 
conditioning is 

PX»*''*|x=x<,Ax = P*^*(Xc). 

For the one-dimensional problem used for illustration 
here, it is easy to compute Pv for a given vector v by 




10 20 30 40 50 



Cumulative timestepper evaluations 

FIG. 4: (Color online) Convergence of Newton-GMRES based 
on the coarse timestepper for a stable steady state (point 2 in 
Fig.EJ. The y-axis is the norm of the (coarse) residual. The 
dashed line shows the estimated magnitude of the stochastic 
simulator fluctuations rescaled to account for replica averag- 
ing (which we use to estimate the fluctuations in the evalua- 
tion of the coarse timestepper). The relative error of the cor- 
responding solution guess at each Newton step is shown in the 
insets (solid (dashed) line for U{V) with y-axis ranging from 
-0.5 to 0.5 except the first one) along the convergence curve. 
The converged coarse steady solution profiles are shown at the 
upper right. Parameters used: m = 40 boxes in the domain 
[0, 0.2], T = 0.05, Lu = 2000 with 150 copies. 

solving P^^y = v through direct linear algebra (e.g. 
Gauss elimination). In general, however, it is worth not- 
ing that P~^y = V can be solved for y with GMRES, 
through repeated calls to the (deterministic) timestep- 
per of the (not-so- accurate) deterministic PDEs (Eqs. ^ 
and 0). 



V. RESULTS AND DISCUSSION 

For our coarse timestepper (based on SSA simulation), 
we discretize the one-dimensional domain [0, L] with 
L — 0.2 into m — AO equally spaced boxes. We choose 
Lu = 2000, which means that the unit density in each box 
corresponds to 2000 molecules. Each evaluation of the 
coarse timestepper corresponds to evolving 150 replicas 
of the SSA simulator for time T = 0.05; the average of 
all replicas is reported. 

We have assumed that there exist some "underlying 
PDEs" that describe the evolution of the (statistics of 
the) SSA simulator averaged over an infinite number 
of copies. In our computations we estimate the coarse 
timestepper of these "underlying PDEs" by averaging 
over several (here 150) copies; even though averaging 
reduces their variance, fluctuations will always remain 
when a finite number of copies is used. Note that we 
should not, in general, expect these "underlying PDEs" 
to be the same as the mean field PDE system described 
in Section ITTll (Eqs. © and lO), which corresponds to 
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FIG. 5: (Color online) Convergence of Newton-GMRES 
based on unpreconditioned (equation-free) and precondi- 
tioned (equation-assisted) coarse timestepper for an unstable 
coarse steady state (point 4 in Fig. |3J . The dashed line shows 
the estimated magnitude of the stochastic simulator fluctua- 
tions, reseated to account for replica averaging (which we use 
to estimate the fluctuations in the evaluation of the coarse 
timestepper). The upper right insets show the converged 
coarse steady solution profiles. Parameters used: m = 40 
boxes in the domain [0,0.2], T = 0.05, uj = 2000 with 150 
copies. 
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FIG. 6: (Color online) Comparison of the leading coarse eigen- 
values of Jacobian matrices based on deterministic timestep- 
per X'^''*'^*, the coarse timestepper D^^*, and the coarse 
timestepper after preconditioning PX'^'^*, which are all eval- 
uated at the computed unstable coarse steady state, (a) 
Leading 40 (smallest magnitude) eigenvalues of X)^"^'*. The 
eigenvalues have already started clustering at 1. (b) Lead- 
ing 40 eigenvalues of T)^"* . This (partial) spectrum is very 
similar to the spectrum in (a). The appearance of complex 
conjugate eigenvalue pairs close to 1 and eigenvalues larger 
than 1 as in (b) (and possibly also in (c)) is probably caused 
by the relatively large fluctuations in the evaluation of the 
coarse timestepper (while the differences between the eigen- 
values within this cluster are relatively small), (c) Leading 15 
eigenvalues of PX'^''*. Most of the eigenvalues are clustered 
close to 1. 



the limit of infinite numbers of molecules. We do, how- 
ever, know that if the parameter lo increases to infinity, 
the "underlying PDEs" do converge to Eqs. ijHl and Q. 

We start by using matrix-free Newton-GMRES to com- 
pute representative spatially nonuniform (both stable 
and unstable) coarse steady states through the SSA sim- 
ulator. The Newton-GMRES fixed point solver used is 
adapted from the MATLAB code nsot?^ with two mod- 
ifications: (a) constant relative tolerance for each GM- 
RES solve; and (b) since our preconditioner changes ev- 
ery time we update the current solution guess, an ad- 
ditional function for constructing the updated precon- 
ditioner was included as an additional input parameter. 
The convergence of Newton-GMRES to a coarse, spa- 
tially nonuniform, stable steady state (point 2 in Fig. O 
is shown in Fig. 21 The magnitude of the fluctuations 
introduced by the stochastic simulator is estimated at 
this coarse steady state as follows: from long-term SSA 
simulations, after the stationary state has been reached, 
we estimate the standard deviation ai for each Ui and 
Vi, i = 1...TO; our estimate of the magnitude of the 
fluctuations is the Euclidean norm of this vector of stan- 
dard deviations. Averaging over n replicas should scale 
this estimate by a factor of ^/n; the resulting estimate 
is marked by a dashed line in Figs. 01 and This num- 
ber also provides an estimate of a reasonable expected 
residual norm upon convergence of the Newton-GMRES. 

Because of the presence of fluctuations in the evalua- 
tion of the coarse timestepper, an inexact matrix-vector 
product is computed at each GMRES iteration. The 
convergence of GMRES in the presence of noise is the 
focus of extensive study in the current literature (see 
e.g. RefsiS2iSii^). These references include discussions 
of bounds of the attainable accuracy of the computed so- 
lution and possible relaxation strategies in the presence 
of noise. In the context of Newton-GMRES, the right 
hand side of the linear equation we want to solve at each 
Newton step (i.e. — VE'(xc)) is also computed with fluctu- 
ations. 

A representative unstable coarse steady state (point 
4 in the bifurcation diagram) is also computed with 
Newton-GMRES. In this case, however, we also im- 
plemented the equation-assisted preconditioning of the 
coarse GMRES, as discussed in Section HV1 using the in- 
verse of the corresponding Jacobian computed from the 
known deterministic approximate PDEs. The conver- 
gence of Newton-GMRES before and after precondition- 
ing (that is equation-free and equation-assisted, respec- 
tively) is compared in Fig. [3 The leading parts of the 
eigenvalue spectra of I?^''^*, IJiif^* (the first 40 eigen- 
values) and PX>'S'** (the first 15 eigenvalues) evaluated 
at the coarse steady state are computed using the coarse 
timestepper and the iterative eigenvalue solver ARPAGK 
(implemented in MATLAB as function eigs) and shown 
in Fig. |S1 A quick inspection of the numerically com- 
puted leading spectra shows that, after preconditioning, 
the eigenvalues of the preconditioned Jacobian based on 
the coarse timestepper were indeed more clustered close 



7 



to 1, consistent with the reduction in GMRES iterations 
observed in Fig. |31 

Note that the first vector in the Krylov subspace con- 
structed for GMRES (vi) is generally obtained by setting 
the initial solution guess for the linear system to zero. 
This implies that vi is different for the preconditioned 
and unprcconditioned GMRES (b and Pb for the linear 
system Eq. (jlOf) . respectively, where P is the precondi- 
tioner). Since all the subsequent vectors in the Krylov 
subspace are built upon the first ones, this may also lead 
to a difference in the number of GMRES iterations to 
convergence. 

Figure [S] shows the cumulative calls to the coarse 
timestepper needed to reduce the nonlinear residual, 
which is similarly defined for both the preconditioned 
and the unprcconditioned case. The results indicate that 
the preconditioner is effective in reducing the nonlinear 
residual, and efficient in terms of overall timestepper eval- 
uations. The preconditioned and unprcconditioned lin- 
ear residuals, on the other hand, are measured in dif- 
ferent norms (||P(b — AAxfc)||2 and ||b — AAxfc||2); we 
impose the same relative tolerance for convergence. For 
non-noisy problems and very tight relative termination 
tolerances, the results of the two types of linear solve 
at the end of the first Newton step would be practically 
the same; with larger termination tolerances, given the 
presence of noise, this is clearly not the case. When the 
initial guess is far away from the true solution (at the 
first Newton step) the initial tolerance for GMRES can 
be set relatively high, to avoid "oversolving" the linear 
equation at the early stages of convergence. 



VI. CONCLUSION 

The purpose of this paper is to illustrate a simple idea: 
that coarse-grained, macroscopic equations can be used 
to assist detailed, fine scale stochastic simulations even 
when they are not really accurate. This is accomplished 
by using certain features of such closed-form macroscopic 
equations (such as their discretized linearizations) as pre- 
conditioners in equation-free iterative linear algebra com- 
putations. This is then an "equation-assisted" approach: 
we compute with a coarse timestepper based on the fine 
scale model, but accelerate the convergence of these com- 
putations using "the best available" continuum determin- 
istic model. 

In this paper we illustrated the concept using a coarse 
timestepper based on a spatially distributed SSA reaction- 
diffusion implementation of the Schnakenberg kinetic 
scheme, and preconditioning with the corresponding Ja- 
cobian derived from the mean-field PDEs. This al- 
lowed us to accelerate the equation-free computation of 
both stable and unstable spatially structured reaction- 
diffusion steady states. The approach can be used as 
a computational "wrapper" around different types of in- 
ner stochastic simulators. The inner simulator was based 
on spatially discretized SSA; the approach could also be 



wrapped around "already accelerated" SSA schemes (e.g. 
those exploiting separation of time scales22i2i^) . The 
approach could also be wrapped around non-SSA, lattice 
gas spatially distributed kinetic Monte Carlo simulators, 
or around simulators based on "the best available" an- 
alytically coarse-grained models of kinetic Monte Carlo 
processes (e.g. Refs i^^i^^ ). It can also be wrapped around 
different (non-kMC) types of fine scale or hybrid models 
such as Lattice-Boltzmann inner simulatorgifl (with den- 
sity PDE preconditioning), or around molecular, Brow- 
nian or dissipative particle dynamics simulators of con- 
densed matter problems, with the preconditioning com- 
ing from traditional continuum closures (elasticity the- 
ory, non-Newtonian rheology). Beyond steady state com- 
putations, such preconditioning might also be helpful in 
other coarse-grained computations involving matrix-free 
iterative linear algebra, such as implicit coarse integra- 
tion schemes. 



VII. ACKNOWLEDGEMENT 

This work was partially supported by the U.S. Depart- 
ment of Energy (IGK, LQ through PPPL) and DARPA, 
by the NSF (DMS-0404537, CTK) and by the Biotech- 
nology and Biological Sciences Research Council, Univer- 
sity of Oxford and Linacre College, Oxford (RE). It is a 
pleasure for the authors to acknowledge helpful sugges- 
tions by Giovanni Samaey and Wim Vanroose during the 
preparation of this manuscript. 



APPENDIX: THE BIFURCATION DIAGRAM 

The bifurcation diagram with respect to the diffusion 
coefficient c?2 of species V is computed using both the 
steady state and the deterministic timestepper approach 
with identical results and is plotted in Fig. El 

We used the following parameters in our computations 
with the deterministic timestepper from the discretized 
ODE system: domain length L = 0.2, number of nodes 
TO = 40, time reporting horizon T — 1. 

The stability of the steady state solutions is identified 
by checking the leading eigenvalues (A^) of the Jacobian 
matrix of the linearized ODE system evaluated at the 
steady states; we confirmed that the eigenvalues of the 
linearization of our timestepper at steady state {^i) in- 
deed satisfy Xi = hi fit /T. The matrix-free eigenvalue 
solver ARPACK is used to compute the leading eigenval- 
ues for both the deterministic and the coarse timsteppers. 

The first two Fourier coefficients of the steady state 
solution, a and b (for sin(a::) and cos(x) respectively), 
are plotted versus the bifurcation parameter d2. The 
(partial) bifurcation diagram consists of four different 
branches of solutions (plotted in different colors) . 

The steady states computed in this bifurcation dia- 
gram show at most one peak due to the relatively short 
domain length. Steady states with n peaks (as shown in 
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Fig. ^ can be easily obtained by using n copies of the 
one-peak solution as building blocks; yet the stabilities 



of the multi-peak and one-peak steady states are not the 
same (see e.g. Refi^). 
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